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Abstract 

We propose an estimator denned in real space for the reconstruction of the weak lensing potential 
due to the intervening large scale structure from high resolution maps of the cosmic microwave 
background. This estimator was motivated as an alternative to the quadratic estimator in har- 
monic space to surpass the difficulties of the analysis of maps containing galactic cuts and point 
source excisions. Using maps synthesised by pixel remapping, we implement the estimator for two 
experiments, namely one in the absence and one in the presence of detector noise, and compare the 
f^ reconstruction of the convergence field with that obtained with the quadratic estimator defined in 

harmonic space. We find good agreement between the input and the reconstructed power spectra 
using the proposed real space estimator. We discuss interesting features of the real space estimator 
and future extensions of this work. 
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I. INTRODUCTION 

Observations of the cosmic microwave background (CMB) are a powerful probe of the 
physics of the early universe and a robust discriminant of cosmological models. Since the 
CMB photons propagate more or less freely after recombination, what we observe on large 
scales are mainly the perturbations on the last-scattering surface. The CMB contains also 
an imprint of the photon quadrupole which developed during the recombination epoch and 
generated linear polarisation through the mechanism of Thomson scattering. 

However, the CMB should also contain signatures of processes that occurred between 
the last-scattering surface and the present. These processes will be manifest as non-linear 
signals which become important on small scales. The growing improvement in resolution 
and sensitivity of CMB experiments allows us to resolve ever so smaller scales and motivates 
the development of estimators and robust implementations for the extraction of the different 
signals. Among such processes is weak gravitational lensing, which consists of the deflection 
of CMB photons by the mass of matter clustered along the line of sight, while leaving 
the surface brightness unaffected. This effect distorts the temperature-temperature (TT) 
correlation power spectrum by a few arcminutes but coherently over degree scales jl] , which 
means that lensing becomes an important effect for £ > 3000. Weak lensing also mixes 
and distorts the EE and BB polarisation power spectra [2]. Finally, weak lensing also 
correlates with other secondary effects, such as the integrated Sachs- Wolfe (ISW) effect and 
the Sunyaev-Zel'dovich (SZ) effect, thus inducing non-Gaussianities which are manifest in 
higher-order point correlation functions [3111]. 

The relevance of the CMB lensing reconstruction for cosmology is threefold. First, by 
converting a fraction of the dominant E-mode polarisation due to density perturbations into 
a B-mode polarisation, CMB lensing introduces a contaminant in the measurement of the 
primordial B-mode polarisation and consequently of the inflation energy scale [5]. Second, 
CMB lensing is also a powerful cosmological probe of the matter distribution integrated 
from the last scattering surface at redshift z ~ 1000 to the present time. Third, by changing 
the power spectra of the perturbations and inducing non-Gaussianities, CMB lensing affects 
cosmological parameter estimation. Delensing the CMB will thus allow one to recover the 
primordial B-mode, probe the full-sky large-scale structure distribution with a maximum 
efficiency at z ~ 3 [0] , and obtain unbiased parameter constraints from CMB measurements 

Effort has been devoted to developing an optimal reconstruction of the lensing potential 
in harmonic space, which implicitly assumes full-sky coverage without galactic cuts, bad 
pixels due to excision of point sources or nonuniform weighting for uneven sky coverage 
[8]. In this idealized context, it has been shown how to construct an optimal quadratic 
estimator for a reconstruction based on the temperature anisotropy alone [9] (see [10] for an 
extension of the formalism to include the polarisation). The improvement gained from using 
a maximum likelihood, and thus more optimal, estimator for the reconstruction from the 



temperature is marginal [TT] (although this might not necessarily be the case for experiments 
sensitive to the B-mode polarisation [12]) . This is because the distortion due to lensing is 
small compared to the intrinsic cosmic variance or to the experimental noise. 

Our interest lies in considering a slightly less optimal estimator modified to have a finite 
range in real space. The estimator is based on a convolution of the square of the lensed 
temperature map with a local kernel, as has been suggested previously [13]. An estimator 
with kernel in real space acts locally, as opposed to the conventional estimator with kernel in 
harmonic space which acts over the entire celestial sphere even though it has finite support 
in harmonic space. The extent in real space will be the inverse of that in harmonic space, 
thus finite and determined by the scale of the root mean square of the lensing convergence, 
with the contribution from sources beyond the kernel being negligible. An estimator with 
these properties can easily be adapted to include cuts, excisions of pixels and nonuniform sky 
coverage for a local extraction of the lensing convergence. In particular, a nonuniform sky 
coverage can be accounted for by a corresponding weighting in the real space convolution. 

The lensing distortion of the CMB anisotropy on the surface of last scattering can be 
described in three ways: by the lensing potential if), by the deflection vector ex = Vip or by 
the convergence tensor k, = — VV^/2. The descriptions of if) and ex suffer from an ambiguity 
upon translation, since a patch of the sky and its translation have the same likelihood on 
account of isotropy. In contrast, the description of the convergence, which is a gradient of 
the deflection vector field, is locally well defined. Here we implement the real space estimator 
to reconstruct the convergence and compare the reconstruction to that obtained with the 
harmonic space estimator. 

The convergence tensor can be decomposed as a sum of an isotropic (diagonal) tensor 
and an anisotropic (traceless) tensor as follows 

K=( K0 + K+ KX V (1) 

For a small patch of sky, we can use the flat-sky approximation and define a Cartesian 
coordinate system with basis vectors {e x (r), e y (r)} and metric g = diag(l, 1). In particular, 
the convergence kq = —(0% + dy)ip/2 magnifies a feature on the last-scattering surface, the 
shear k + = — (<9^ — <9^)-?/>/2 stretches it along the x-axis while compressing it along the y-axis, 
and analogously the shear k x = —d x d y ip stretches it along the y = x axis and compresses 
it along the y = — x axis. Here we focus on the dilation effect captured by the isotropic 
component Kq, and will treat the shear components in a forthcoming study. 

The convergence can be further decomposed into two contributions, namely the conver- 
gence of the unlensed temperature anisotropies, to a high degree isotropic, and the con- 
vergence of the lensing field, by nature random. In order to remove the zero mode of the 
unlensed CMB, we average the convergence on each local patch over the sky region surveyed 
and subtract this average value from the estimated convergence map. Higher multipoles 
of the unlensed CMB will still be present in the estimated convergence map, which give 



rise to a noise bias. This bias can be removed by cross-correlating the reconstructed map 
with maps of the large scale structure derived from independent methods, e.g. a galaxy 
redshift survey. It can also be removed when estimating the convergence power spectrum 
by combining, in harmonic space, different triangle configurations with a common side. The 
formulation of an estimator for the lensing convergence given the lensed image still requires 
some prior knowledge of the statistics of the underlying temperature field, in particular the 
power spectrum of the unlensed sky which is assumed to be Gaussian. 

This is the first of a series of studies where we explore the properties of the proposed 
estimator and test its implementation. The manuscript is organized as follows. In Section [TT] 
we describe the weak lensing CMB of photons in the Born approximation. In Section |IH| 
we derive the proposed estimator defined in real space, whose kernel in harmonic space is 



described in Appendix [A} In Section IV we proceed to reconstruct the convergence power 
spectrum by implementing the real space estimator on lensed maps synthesised as detailed 
in AppendixlBJ For comparison, we also implement the harmonic space estimator as detailed 



in Appendix O We discuss the results obtained from the two estimators in Section V and 



conclude in Section VI by discussing further developments of the real space estimator in 
follow-up studies. 

II. WEAK LENSING OF THE CMB 

Weak lensing of the CMB consists of the deflection of CMB photons from the original 
propagation direction 6 on the last scattering surface (the source plane) to an observed 
direction on the sky today (the image plane). This deflection amounts to remapping the 
unlensed temperature anisotropies T to the lensed ones T according to 

f (0) = T{6). (2) 

The deflection angle a = — is given by a = Vip, where the lensing potential ip is given 
by the gravitational potential projected on the image plane and integrated along the line of 
sight, and V is the covariant derivative on the image plane. The lensing potential measures 
the cumulative effect of deflectors along the line of sight which lens the perturbations at the 
last scattering surface and thus lensed are imprinted on the CMB. Consider the line element 
of a homogeneous and isotropic (l+3)-dimensional spacetime subject to scalar perturbations 
expanded up to linear order 

ds 2 = a 2 ( V ) [-(1 - 2^)dr, 2 + (1 + 2*) (dr 2 + f 2 (r) dtf k )] , (3) 

where the functional form of fk depends on the spatial curvature k of the two-dimensional 
surface Qk as follows 

for k = +1, 
fk(r) ■= { r : for k = 0, (4) 

for k = —1. 




Integrating the equation for a null geodesic in an arbitrarily curved spacetime, we find for 
the lensing potential ip generated by the gravitational potential \l/ that [6] 

^-'T^wm^' (5) 

where r is the conformal distance to the lens plane, r^s is the conformal distance to the 
last scattering surface and 6 is the two-dimensional position vector on the image plane. 
This expression was derived by integrating the photon geodesies from the source along the 
unperturbed path (the Born approximation). The gravitational lensing by ip deflects the 
photons on the image plane by V"0, thus remapping the temperature anisotropics according 
to 

f (0) = T(0 + W) = T(0) + V^ ■ VT(0) + 0[(V^) 2 ]. (6) 

Since the Born approximation assumes that Vip is constant between and 0, it is only valid 
for deflections small compared to the scale of the lensing perturbations, i.e. deflections up 
to a few arcminutes. 

Considering the full sky described as a two-sphere centred at the observer, then the 
temperature anisotropies can be expanded in spherical harmonics as follows 

oo e 
T W = J2 J2 a *™ Y e™(ttk), (7) 

1=0 m=-t 

where denotes the position on the two-sphere. In the flat-sky approximation, the CMB 
temperature anisotropy is instead expanded in terms of plane waves as follows 

1 oo 

T{0) = —j= V a e exp[i£ ■ 6] (8) 

\/47T ^"^ 

where 6 = (8 x ,6 y ). To obtain the same asymptotic density of states as on the full-sky 
description, we compactify the sphere to a square having sides of length L = y^n so that 
the discretized variable 0/L is replaced by the continuous wave number £/(2tt), yielding the 
Fourier mode expansion 

This is a valid approximation for sufficiently small angular scales where we can replace the 
expansion in spherical harmonics by an expansion in plane waves. At the largest angular 
scales, however, the curvature and connectivity properties of the celestial sphere are no 
longer negligible and the approximation is no longer valid. Likewise expanding the lensing 
potential in plane waves so that (Vijj)(£) = i£ip(£), we can write the linear order correction 



of the lensed temperature anisotropy in harmonic space as the convolution of the lensing 
potential with the unlensed temperature anisotropy 

no =w>- J *"■(*-*) mm-*) +<&(*>]. do. 

The lensing effects are apparent in the temperature power spectrum. Defining the tem- 
perature power spectrum by (T(£) T(£')) = (2tt) 2 S(£ + £') Ce, we find that the power 
spectrum of the lensed temperature anisotropy is related to that of the unlensed as follows 

f{£') f{£ -£')) = (2tt) 2 S(£) C e + (2tt) 2 [£ ■£' C e , +£■(£- £') C\ t _ t >\] ^(£). (11) 

Since the convolution is made of gradient terms, it results in the transfer of power from large 
to small scales, with the random deflections simultaneously smearing out the sharp features 
of the unlensed spectrum and thus leading to a suppression of the acoustic oscillations 



PP . The first term in Eqn. (11) describes a statistically isotropic ensemble where modes 
with different £ are uncorrelated. Anisotropies, and in particular those generated by the 
lensing potential, introduce correlations among different £ modes. What we will be using to 
reconstruct the lensing field are these off-diagonal correlations. 



III. ESTIMATOR OF THE LENSING CONVERGENCE IN REAL SPACE 

A. Derivation of the real space estimator 

In harmonic space the estimator of the lensing potential ip(£) is expressed as the convo- 
lution of the square of the lensed map T(£) by a weight function Q 1 ^ (£,£') as follows 

m = J-0p f ( £ ') f ( £ - £ ') &&*)■ (12) 

The optimal estimator is found for the weight function which minimizes the leading order 
variance, which is given by 

hl „ r l£-£'C e , +£■(£- £')Cu £ i\ , s 

Q*(l,e)=Nt - r ^ LJ^H ( 13 ) 



where N( is the detector noise power spectrum given in Eqn. ( A10 ) and He, is the variance of 



the estimator given in Eqn. (A13) (see Appendix [A] for a derivation). The relation between 



the convergence and the lensing potential in harmonic space translates into the same relation 
for the corresponding estimators. Hence the weight function for the convergence estimator 
will be 

Q(£,£')=£ 2 Q^{£,£'). (14) 



Naively we would write the corresponding estimator in real space as 



= fd 2 e'f(e') fd 2 0"f(0") Q(0,0',e"), (15) 

where we define the corresponding weight function in real space by 

Q(0, 0', 0") = J -0- exp[i£ ^]j-0y 2 exp[-^' • 0'] exp[-i(i - (!) ■ 0"] Q(£, £'). (16) 

In general, Q(0, 0' , 0") would be a function of six variables, namely the three lengths and 
the three corresponding angles of the three vectors 0, 0' and 0". However, the functional 
dependence of Q(0, 0' , 0") is derived from Q(£, £') which depends on three variables, namely 
the lengths of £ and £', and the angle between them £/ = <p£ — (fig, due to the fact that the 
vectors in harmonic space form a triangle. This reduces Q(0,0' ,0") to a function of three 
variables, namely the length of two of the vectors and 0' , and the angle between them 
£,o = <Pe ~ 00' • We can expand the kernel in terms of eigenfunctions that factorize the radial 
and the angular dependence as follows 

+oo 

Q(£,£')= Y, exp[im&]Qm(40, (17) 



m=— oo 



so that 



QM?) = ^J d & e M-i™&] Q(£,£'). (18) 

Let (pi and <\>g be the angles that £ and £' respectively make with the £ x -axis, and likewise 
(pe and (f) e i the angles that and 0' respectively make with the 8 x -a.xis. By construction 

£' .0' = M COB[fo - fa], (19) 

(£ - £') ■ 0" = W" cos[<^ - 00"] - l'e" cos{4>t! - <po"] • (20) 



Without loss of generality, let us also consider = 0. Then Eqn. (16) becomes 

+0O 

Q(0',0",6O = J2 exp\im4o]Q m (0',0",Zo), (21) 

m=— oo 

where 



-I POO /'OO 

QUO', 0", 60 = jtt^ / dll d£' £' J m (i0") Jmit'ie' - 6" cosfo,])) Q m {t, £') (22) 

I/ 71 " J Jo Jo 

and £0 = / — » is the angle between 0' and 0" . 1 The eigenfunctions Q m (9', 9", £ ) depend 
explicitly on £g which means that the factorization of the angular and radial dependence of 



1 Here we used the Jacobi- Anger expansion of the Bessel functions [TJ] 

oo oo 

exp[i£9 cos[</>]] = ^ i n J n (iO) exp[in<f>] = J (£6) + 2 ^ i n J n (£8) cos[n</>] (23) 



n— — oo n— 1 

7 



the kernel in harmonic space failed to translate into the same factorization of the kernel in 



real space. The computation of the estimator in Eqn. (15) requires a quadruple numerical 
integral 

/„ +00 

d 2 0' d 2 0" J2 ex P [tm^]f(0')f(6")Q m (9',9"^e) (25) 



m=— 00 



which involves a number of operations proportional to the number of components of Q m , 
hence cubic in the size of Q m . In order to minimize the computational costs, we would like 
instead to write the real space estimator in such a form that the radial eigenfunctions in 
real space depend only on the length of the vectors and not on the angle between them. 

To achieve this we change to the variables £ + ,£- such that £ = £ + , £' = (£ + + £-)/2. 
In the coordinates (£ + ,£_) the kernel becomes squeezed in the £_ direction as £ + increases. 



Then Eqn. (12) yields 



^+) = l^f (^y^) f (^y^) w ^ £ -^ < 26 ) 



d 2 £. 

(27T) 



where 



W(£ + ,£-) = Q(£,£') (27) 

upon change of variables. The corresponding estimator in real space becomes 

k o (0) = I d 2 0+ I d 2 0_ f{0 + + ) f{6+ - ) W{0, 0+, 6 _), (28) 

where we define the corresponding weight function in real space by 

W(e,0+,0J) = ||^ exp[z£ + -e]J^ exp[-i(£ + -0 + + £„-6-)]W(£ + ,£-). (29) 

Again we expand the kernel W(£+, £-) in terms of eigenfunctions that factorize the radial 
and the angular dependence such that 

W m (i + ,iJ) = ±- J d X , expi-imxi] W(£ + ,£„) (30) 

and xt is the angle between £ + and £_. Let 4>i + and 0£_ be the angles that £ + and £_ 
respectively make with the £ x -axis, and likewise <pe + and <p9_ the angles that + and 6 
respectively make with the ^-axis. By construction 

£ + 6 + +£ 6 = t + + cos[<l> l+ - <p e+ ] + £_#_cos[0 £ _ - <f> 9 _\. (31) 



and the orthogonality condition 

J d<j) exp[i£6cos[(j>]]exp[in'<p] = l dcj> ]T i n J n {£9) exp[i{n + ri)<j>] = {2n)i n ' 5 nn , J n {l0). (24) 



71 — — 00 



For the case that = 0, Eqn. (29) becomes 



+00 



W(6 + ,6_,xo) = J2 eMimxe}W m (9 + ,9. 



(32) 



m=— 00 



where 



WmW+,6- 



.1. J m {i+B + ) J m (i-e-) W m (£+,£-) (33) 



(2^) 2 Jo T jo 
and xe = 4>e + — 4>o_ is the angle between + and 0_ . We have now achieved the factorization 



of the angular and the radial dependence of the kernel in real space, so that Eqn. (28) can 
be written as 



«o< 



+00 » „ 

= 0) = J2 / d2 °+ / d2 °- eM* m Xe] f(8+ + 0_) f (0+ - 9 J) W m (9 + , B-). (34) 



m=— 00 



In contrast with Eqn. (25), this integral involves a number of operations proportional to the 



square of the size of W m . 

The convergence map thus estimated contains two contributions: the convergence derived 
from the lensing potential and a gaussian noise derived from the unlensed CMB, k = 
ftolv + Ko\ij>=o (as we will see in Section [Vl the (white) detector noise averages out for our 
pixel based estimator). Taking the average of the convergence map (with a fixed lensing 
potential) over many realisations of the CMB results in the gaussian CMB noise averaging 
out so that our estimator is unbiased (k ) = («olv) • However, since in practice we only have 
access to a single realisation of the sky, our reconstructed convergence map will be noisy. 



In the absence of lensing, the convergence estimator defined in Eqn. (34) will have a non- 



vanishing average value. We therefore subtract the average value of the map reconstructed 
in the absence of lensing so that our convergence estimator becomes 



«o(0 = 0) 



d 2 6 + / d 2 6 ^ exp[imxe] 



+00 

E 

m=— 00 ' 

x [t(0 + + 6L)T(6L 



6 ) - (f(0 + + 0_)f(0 + - 0_)\ W m (0 + , #-), (35) 



where the average is taken over many realisations of the unlensed CMB. 

We are interested in estimating the power spectrum C^° K ° from the reconstructed conver- 
gence map. In the absence of lensing there is a non-vanishing contribution to the estimator 
from the collapsed four point function (k |^=o k oU=o) °^ ^ ne unlensed temperature field, 
which biases the power spectrum of the estimated convergence. This bias is precisely the 
variance of the gaussian CMB noise in the reconstructed convergence map [15J. To remove 
this bias we generate a set of unlensed temperature maps (based on the unlensed tem- 
perature power spectrum) and compute the convergence power spectrum for each of these 
maps. The average of these power spectra is then subtracted from the power spectrum of 



the estimated convergence to yield an unbiased estimate of the convergence power spectrum 

(k Q «q) = (ftolv> k oU) • 



In this formulation, the kernel in Eqn. (35) is essentially the Bessel transform of the 



optimal weight function in harmonic space. Contrary to the formulation in harmonic space, 
where the temperature map is transformed before being weighted by the kernel in harmonic 
space, here it is the kernel that is transformed before acting on the temperature map to 
yield the lensing potential. The kernel weights the contribution of each pair of points on the 
temperature map to the weak lensing. Since weak lensing is manifested essentially at very 
small scales in the temperature map, all the information relevant to lensing reconstruction 
lies on angular scales close to the resolution scale of the temperature map, which is set 
by the beam size for both the harmonic and real space estimators. In the case of the real 
space estimator there is an additional effect that limits the reconstruction on small scales, 
which arises from the fact that we are averaging the squared temperature map over the 
spatial extent of the kernel. Since the kernel changes sign, the averaging over quadratic 
combinations of closely separated pixels results in a loss of recovered power on angular 
scales smaller than the spatial extent of the kernel. This effect has been studied in more 
detail in Ref . [16] . The angular resolution of the experiment limits the spatial extent of the 
kernel so that an experiment with higher angular resolution or higher sensitivity will have a 
more local kernel. 

IV. RECONSTRUCTION OF THE CONVERGENCE MAP 

We proceed to describe the implementation of this estimator to reconstruct the conver- 
gence field kq- First we synthesise lensed CMB maps as described in Appendix [Bj We then 
implement the real space estimator to the maps by applying the procedure described be- 
low, and produce the maps of the reconstructed convergence field. For comparison of the 
performance of the real space estimator, we also implement the harmonic space estimator 
by applying the procedure described in Appendix [Cj In this study we will be using two 
experiments inspired by the specifications of the PLANCK experiment for the v = 143 GHz 
channel [17J and only differ in the detector noise. The experiments, denoted by Designer 
and Planck, both have 6f w hm = 7.8 arcmin for the beam full-width at half maximum, 

1 I s } 

and a pix = and a pix = 6.8 / / fiK/r&d for the white noise amplitude per beam width 
respectively, where f s k y is the fraction area of the sky covered by the patch [TBI US]- We 
chose the experiments to have the same beam size in order not to entangle the effect of the 
noise with that of the beam. In future work we will study experiments with higher angular 
resolution such as ACT |20| and SPT 
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Figure 1: Contour plots of Q m (£,£'). The levels denote the orders of magnitude as fractions 
of the maximum value, delimited by {1, 1CP 1 , 1CP 2 , 10 -3 , 10 -4 } from the lightest to the darkest 
shade. Here m max = 16, £ max = 4000 for PLANCK noise at v = 143 GHz. 

A. Generate the kernel 

The independence of the kernel eigenfunctions on the angular coordinate in the coor- 
dinates (£+,£.-) renders this parametrization more attractive to implement. We start by 
defining the function W (£+,£-) given by Eqn. (27). We then compute W m (£ + ,£-) using 
Eqn. (30). The corresponding inverse Fourier transform W m (0 + ,0J) is computed using 



Eqn. (33), with the final expression for the kernel following from Eqn. (32) 



Two parameters that we need to define for the computation of the kernel that affect the 
reconstruction of the lensing convergence are the kernel extent in I space, £ m ax, and the 
number of eigenfunctions at which we will truncate the series in m, m max . Since the kernel 
is computed from the lensed map, the value of £ max should be determined by the range in £ 
space that the experiment can probe before the beam and noise cut off the signal. We thus 
set £ max equal to the higher £ moment that we use to synthesise the maps, i.e. £ max = 4000. 
This is also consistent with the value read off of Q m (£, £') which we plot in Fig. d] 2 We have 



The small scale structure seen in Fig. [I] is likely numerical noise due to the discretization. Various tests 
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Figure 2: Contour plots of W m (@+, &-)• The levels denote the orders of magnitude as fractions 
of the maximum value, delimited by {1, 10 , 1CP 2 , 10 -3 , 10 -4 } from the lightest to the darkest 



shade. Here m r 



4000 for no detector noise. 



tried using higher values for £ max and found that neither the kernel nor the reconstructed 
power spectra changed significantly due to the fact that most of the lensing information is 
accessed around the beam scale. In the absence of detector noise, the kernel is displaced 
to higher values of £ and £', up to the chosen £ ma x- To fix the value of m max we compared 
the relative amplitude of each eigenfunction W m (0 + , 0J) and discarded those that were less 
than 1% of the dominant m = eigenfunction. This corresponds to m max = 4. 

For each experiment we generate W m (9 + , 6 J) for the selected values of these two param- 
eters. In Figs. [2] and [3] we show the contour plots of four eigenfunctions of the kernel in real 
space for the Designer and the Planck experiment respectively. From the plot of the 
m = eigenfunction (which we have found to dominate the lensing reconstruction) we read 
off the value for the spatial extent of the kernel 9 max , used in the convolution routine to 
compute the convergence map. In harmonic space the kernel peaks at the scale of the beam 



that we tried seem to indicate that the quality of the lensing reconstruction is not significantly affected by 
this small scale structure. Nevertheless, we plan to improve the numerical calculation of the kernel and 
its transformation in future higher precision studies. 
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Figure 3: Contour plots of W m (@+, &-)• The levels denote the orders of magnitude as fractions 
of the maximum value, delimited by {1, 1CP 1 , 1CP 2 , 10 -3 , 10 -4 } from the lightest to the darkest 



shade. Here m r 



4000 for PLANCK noise at v = 143 GHz. 



but has support at smaller £ values. The spatial extent of the kernel is thus larger than the 
beam scale and scales roughly linearly with the beam size. We observe that the kernel is 
compact and extends only up to 0.70° and 0.60° at the 1% level in the absence and in the 
presence of detector noise respectively. In practice, the cut-off in the kernel was tested by 
comparing the power spectrum of the convergence map reconstructed for different values of 
9 max . For each experiment, we verified that increasing the value of 9 max beyond 0.35° did 
not improve the reconstruction of the power spectrum, whereas reducing the value of 9 max 
below 0.35° degraded the reconstruction. We thus set 9 max = 0.35°. 

B. Compute the convergence map 

We compute the estimator by convolving the kernel with the lensed CMB map according 



to Eqn. (34). We implement the convolution using the discrete version of this equation as 
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Figure 4: Relative sizes in the reconstruction of the convergence map. The outer black 
square represents the map and the inner, dashed square centred at the origin represents the kernel. 
The red square around the kernel encompasses the points on the map which contribute to the 
convergence on the point where the kernel is centred. The kernel moves over the map by successive 
translations indicated by the vector Oq. At the new centre, a new red square is defined. This 
procedure is repeated until the red square reaches the boundary of the map. 
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(36) 



The size of the kernel determines the extent of the temperature map over which the convo- 
lution is computed. Since the temperature map is evaluated at the sum and difference of the 
contributing pixels, we need a region of the temperature map that is four times the length 
of the kernel to compute the convolution (see Fig. 111). At an arbitrary point 6 = (9o x ,9o y ) 
the convergence is given by the estimator at the origin translated by the vector 0. 

We then average the reconstructed convergence map over all pixels to extract the zero 



mode, which is then subtracted from the map as discussed in the text following Eqn. (35). 



Finally we subtract the convergence power spectrum of a Monte Carlo realization of one hun- 
dred unlensed CMB maps from the convergence power spectrum of the lensing potential, 
thus removing the gaussian bias from the estimated power spectrum. Our implementation 
of the real space estimator takes a little over one minute in total on a desktop computer, 
including the calculation of the kernel (about 30 sec) and the reconstruction of the con- 
vergence map (about 40 sec). This is about six times the time that the harmonic space 
implementation takes to reconstruct the lensing field. 
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V. RESULTS AND DISCUSSION 

For each experiment we plot the output power spectrum C^ outKout , both before (light grey 
solid line) and after (dark grey solid line) the removal of the bias, against the input power 
spectrum C^ mKin (black solid line). In Figs. [5] and [6] we plot the reconstructed power spectra 
for the Designer and the Planck experiment respectively. We also plot the variance of 
the estimator Vax[k (l)} (black dot-dashed line) as well as Vai[k (£)}/ f s k y (grey dot-dashed 
line) as an auxiliary tool for visualization. The error bars on the recovered power spectra 
are computed from the standard deviation of the reconstructed power spectrum, binned over 
logarithmically spaced intervals in £ space and added in quadrature. The variance of the 
reconstructed power spectrum encompasses both the variance of the optimal estimator and 
the average value of the convergence fluctuation in the map, weighted by the fraction of the 
area of the sky covered by the patch, f sky , as follows 



Var[C7f K0 ] « — — (Vax[«o(£)] + C?° K °) 2 . (37) 

" Jsky 

We observe that the degree of recovery of the input power spectra is not the same throughout 
the interval in £. In particular, at small £ the cosmic variance dominates the error of the 
reconstruction, whereas at large £ the variance of the estimator takes over. The region where 
Q Ki " K ™ > Var[Ko(^)] indicates the interval in £ where the lensing potential can in theory be 
mapped. Accordingly, this is also the region where the reconstructed convergence power 
spectra tracks more closely the input power spectrum. Also visible are the fluctuations 
in the recovered power spectrum which arise from the subtraction of the noise bias term 
estimated from a finite number of realizations of the unlensed CMB map. 

For the Planck experiment, using the real space implementation we achieve a reasonable 
reconstruction of the input power spectrum in the interval where the input signal is larger 
than the variance of the estimator, although on smaller scales, £ > 600, there is a decrease 
in the recovered power. Comparing the reconstructed power spectra for the Designer and 
Planck experiments we observe that the real space implementation is fairly insensitive to 
the detector noise and returns a reconstruction consistent with the input power spectrum. 
We can understand these results as follows. 

As attested by the power spectra of the maps reconstructed with the real space estimator, 
there is a loss of power at small scales beginning around the angular scale corresponding to 
the size of the kernel. This loss of power is a consequence of averaging modes smaller than 
the finite extent of the kernel and hence is an intrinsic constraint of the real space imple- 
mentation. This intrinsic constraint, studied in more detail in Ref. [H], can be overcome by 
shrinking the extent of the kernel in real space, which can be achieved with a smaller beam 
and detector noise, thereby moving the support for the kernel to larger £. The £ interval over 
which the power spectrum can be recovered by the real space estimator is determined by the 
characteristic lengths of the map and the kernel as follows. The lower £ limit is determined 
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Figure 5: Mean square of the convergence map reconstructed with the harmonic space 
estimator (left panel) and the real space estimator (right panel). The smooth black curve 
is the input power spectrum. The light grey and dark grey solid lines are the output power spectra 
before and after the removal of the bias respectively. The error bars measure the total standard 



deviation binned over logarithmically spaced intervals in 
and rrimax = 4, £ ma x = 4000 for no detector noise. 



Here 



' fwhm 



7.8', FOVr 



map 



66.6° 



by the size of the map, since it measures the longest wavelength mode that can be enclosed 
within the size of the map. The upper £ limit is determined by the size of the kernel, since 
it measures the smallest wavelength mode that the kernel can probe, below which there is 
a loss of power due to the averaging of modes smaller than the size of the kernel. Within 
this £ range, the real space estimator seems to perform fairly well both in the absence and 
in the presence of detector noise. Note that the £ range for a reasonable reconstruction 
with the harmonic space estimator is wider, being limited at higher £ by the angular scale 
corresponding to the beam size, since this is the scale which constrains the action of the 
kernel in harmonic space. 

Comparing the reconstructed power spectra in the absence and in the presence of detector 
noise, we observe that the real space estimator appears to be insensitive to the experimental 
noise. This is because the estimate of the convergence in each pixel is given by the sum of the 
product of pairs of neighbouring pixels weighted by the kernel. As a result, the noise, being 
independent in each pixel, is averaged out. To substantiate this result, we reconstructed the 
convergence map from a pure white noise input map using the real space estimator and the 
harmonic space estimator. In Fig. [7]we observe that the power spectrum recovered by the real 
space estimator is about seven orders of magnitude smaller than the input power spectrum, 
whereas the power spectrum recovered by the harmonic space estimator is comparable to 
the input power spectrum. 

Despite our attempt to remove the CMB and the detector noise bias from the recon- 
structed convergence power spectrum, there still remains a bias in the harmonic space im- 
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Figure 6: Mean square of the convergence map reconstructed with the harmonic space 
estimator (left panel) and the real space estimator (right panel). The smooth black curve 
is the input power spectrum. The light grey and dark grey solid lines are the output power spectra 
before and after the removal of the bias respectively. The error bars measure the total standard 
deviation binned over logarithmically spaced intervals in t. Here Qj w h m = 7.8', FOV map = 66.6°, 
and m m ax = 4, i ma x = 4000 for PLANCK noise at v = 143 GHz. 



plementation, as observed in Figs. [5] and [6j This bias presumably arises from a coupling 
between unlensed temperature modes due to the finite size of the map, providing an ad- 
ditional source of lensing which we have not treated here. This mode coupling is absent 
in the reconstruction over the full sky [22] • The real space implementation appears to be 
insensitive to this bias. 



VI. SUMMARY AND FUTURE WORK 

In this paper we implemented a new estimator of the convergence field for the extraction 
of the lensing potential from CMB maps. Our interest was in reconstructing the convergence 
field from CMB maps from which contaminants, such as point sources and the SZ effect, 
had previously been removed using multi-frequency information. 

The new estimator acts locally in real space, thus being able to treat the excision of pixels 
and nonuniform sky coverage in a flexible manner. We implemented the estimator on two 
experimental setups, one without and one with detector noise, based on the specifications 
of PLANCK for the v = 143 GHz channel. For comparison of the performance of the new 
estimator, we also implemented the conventional estimator which acts in harmonic space. 
From the two experiments we learned that by increasing the range of the kernel in £ space, 
and consequently reducing its range in 9 space, we improve the reconstruction of the power 
spectrum at small scales. The finite extent of the kernel is a desirable property of the 
proposed estimator, since in theory it allows the reconstruction of the lensing convergence, 
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Figure 7: Mean square of the convergence map reconstructed with the harmonic 
space estimator (light grey line) and the real space estimator (dark grey line) from a 
pure white noise map. The grey solid lines are the output power spectra whereas the smooth 



black curve is the input power spectrum. Here 9f w h m 
£ max = 4000 for PLANCK noise at v = 143 GHz. 



7.8', FOV r 



map 



66.6°, and m r< 



manifested essentially at very small scales, from a small map of the sky as long as the kernel 
probes sufficiently small angular scales. Even though there is a loss of power on scales 
smaller than the finite extent of the kernel, this effect could be studied using simulations 
to determine a form factor that could be applied to correct the loss of power [TB]. We 
also observed that although the real space estimator is limited, in a statistical sense, to 
be as good as the harmonic space estimator, it has the advantage over the harmonic space 
estimator of being less sensitive to white noise. In practice, however, CMB experiments 
have to deal with correlated noise, which will pose as much of a challenge to the real space 
estimator as to the harmonic space estimator, most likely requiring the construction of a 
more carefully designed kernel in harmonic space that downweights the correlated modes 
before transforming to real space. 

On the basis of the results presented here, we envisage two follow-up studies. In the 
first study we intend to optimize the current implementation of the real space estimator to 
use a kernel capable of probing smaller angular scales in the reconstruction of the lensing 
convergence, as we would expect for the ACT and SPT experiments. As we have demon- 
strated, the real space estimator can be applied to small patches of sky without incurring in 
the serious spectral leakage that affects the harmonic space estimator on rectangular, rather 
than torodial, domains in the presence of red spectra. In real space, the local filter acts 
up to a kernel length away from the edge of the map without leakage of power. We also 
intend to study the effect of a mask due to the excision of bad pixels, e.g. arising from the 
removal of bright point sources, that can cause aliasing of power from large scales to small 
scales. Inpainting has been proposed as a means of interpolating across masked regions that 



18 



preserves the statistical properties of the map [22] , so we plan to investigate how sensitive 
the real space estimator is to complex masks with and without the use of such techniques. 

A straightforward extension will be to develop the analogous estimator for the shear 
components of the convergence tensor and work out how to combine the different components 
for the reconstruction of the lensing potential. Since the shear components do not contain 
the £ = mode, the corresponding estimators will be less sensitive to a bias. The dilatation 
and shear components provide complementary information on the lensing potential but are 
not independent, being related by consistency conditions that are algebraic in harmonic 
space. Since these components probe the monopole and quadrupole sectors respectively, 
the noise associated with these components will be uncorrelated and the reconstruction 
can be harmonized, with the noise being reduced via inverse variance weighting of these 
components. This point is discussed in more detail in Ref. [16]. 

In the second study we intend to extend the implementation of the real space estimator 
to CMB polarisation so as to optimize the reconstruction of the lensing potential from the 
PLANCK data. 
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Appendix A: Derivation of the harmonic space estimator 

In this Appendix we derive the optimal estimator for the lensing potential in harmonic 



space. From Eqn. (11) a crude estimator for the lensing potential would be the £ ^ term 



Ml £') = l T(£')f (£-£>) 

This way we would be estimating the anisotropic component of the convergence which is 
due to the lensing potential only, i.e. without the isotropic zero mode. 
We now seek for the optimal combination 

4>(£) = J0^H£,£')G(£,£'), (A2) 

where the optimal weights G (£,£') are given by 

d 2 £ f d 2 £' 1 



G(£,£') 



a 2 



ip(£,£') 



^yj ( 2 -) 2 ^) 
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(A3) 



so as to minimize the variance, and are such that 

d 2 £ f d 2 £' 



(2tt) 2 J (2tt) 



G(£,£') 



(A4) 



so as to make the estimator unbiased (ip(£, £')) = if)(£)- Here a 2 - , is the variance of ip(£, £') 
which for the case that ip(£) is centred at zero becomes 



i>(t,€') 



M^l)-^l)) W<2,4)-^2) 



(2vr) 4 



f(^) f(£i -£[) f*{£' 2 ) f*{£ 2 -£' 2/ 
1 



[€i ■ £\ C e[ + 4 • (£1 - £[) C l£l _ e{l ] [£ 2 • 4 Ce 2 + £ 2 ■ (£ 2 - £' 2 ) C^-^] 



2 6{£! - £ 2 ) 5{£[ - £ 2 ) 



C^C^-^l 



We then find for the optimal estimator that 



2' 



(A5) 



${l)=Mi 



where 



1 



d 2 £' ~, , x ~ , , x 1 £ ■ £' C e + £ ■ (£ - £') C le _ e n 
j-^ T(£') T(£ - £') . } ; l£ £l 



d 2 £' ! [£■£' C fJ +£■{£- £') C\ t _ g \} 2 
(2tt) 2 2 : = 



(A6) 



(A7) 



CV C\t-t'\ 

Note that the calculation of the variance presupposes a prior knowledge of the unlensed 
power spectrum of the temperature anisotropies. Here we use the same TT power spectrum 
that we use for the synthesis of the CMB map. 

The error associated with the reconstruction of the lensing potential with this estimator 
is given by the variance of the estimator. We can identify the variance of the estimator 
with the power spectrum of the estimator as follows & 2 -j (ee ,^ = (V'(^i)^i) i J *(^2,£ 2 )) = 
(27r) 2 8{£\ —£ 2 ) 5(£' 1 —£ 2 ) Var[0(£, £')]. It follows that the variance of the optimal estimator 
is given by 

-i 
= (2n) 2 5(£! - £ 2 ) Var$(£)] 



a 



i>(i) 



2 ei 



d 2 £ 



(27T) 



where 



Var 



d 2 £[ 1 



d 2 £' ![£■£' C e , +£■{£- £') C^,] 



(Ah 



(2tt) 2 2 



Cfi C\i_i'\ 



Aff. 



(A9) 



The variance gives an upper limit to the dispersion about the mean reconstructed map of 
the output maps reconstructed with the estimator from a finite sample of input synthesised 
maps, thus being a measure that can be used to validate the estimator. 
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The expression for the optimal estimator can most straightforwardly be modified to in- 
corporate detector noise and finite beam width by replacing CV by Cy + Ny. The noise 
power spectrum N £ is given by the inverse-sum over the number of channels of the detector 
noise rii(£) of each channel % 



'num_chann 



w «-(S^)- (A10) 

The detector noise is modelled by a gaussian signal on the beam size and includes both the 
white noise amplitude and the beam profile attenuation factor as follows 

th{£) = (e fwhmi a piXi ) 2 ex P [(9 fwhmi ) 2 £(£ + l)/(8\n2)}. (All) 

Here Qf w h mi is the beam full-width at half-maximum of the beam and u P i Xi is the white noise 
amplitude per beam width. 

We can regard N £ as the total variance of a multiple sampling of size the number of chan- 
nels, each sampling i having variance rii(£). Finally, the estimator for the lensing potential 
becomes 

J (2f) 2 '2 [C,, + N e ]lC lt _ ei + N lc _ ei }' 



where 

1 f d 2 £' 1 [£ ■ £' C ff + £ ■ {£ - £') C\ t _ t ,\Y 



Aft J (2tt)2 2 [C e + N if }[C l£ _ £ri + N l£ _ ei ] 



(A13) 



The estimator for the deflection angle ol(0) = Vip(0) and the convergence ^o(^) — 
—'V 2 ip(d)/2 will be given respectively by the gradient and the divergence in real space of 
i>(£). In £ space this amounts simply to 

a(£) = i£ ij(£) (A14) 

«o(£) = f ${l) (A15) 

with variances Var[o;(£)] = £ 2 Var[-0(£)] and Var[K (£)] = £ 4 Vax[i(j(£)] respectively. 

Appendix B: Synthesis of the lensed CMB map 

For the construction of the CMB map and the lensing map we used power spectra gener- 
ated by CAMB [23J, respectively the TT and the ipip power spectra, on scales with £ < 4000. 
Given the nature of the CMB lensing and the form of the data to be analysed, we will be 
working upon a square patch, thus bounded, of the unbounded celestial sphere. For a proper 
treatment of the boundaries of the patch, we must simulate the unboundedness of the sphere 
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by implementing periodic boundary conditions. In order to reproduce periodic boundary 
conditions, we apodized the CMB map by smoothing to zero over a strip around the bound- 
ary with a given width fraction (in this case 0.10). The requirement of periodic boundary 
conditions makes the implementation of the quadratic estimator in harmonic space easier 
but is not necessary for the implementation of the real space estimator. 

We then applied the lensing map to the CMB map to produce the lensed CMB map. 
This consisted of shifting the CMB map by the lensing potential assuming the Born approx- 
imation, as detailed in the Appendix of Ref. [7]. A fast way to construct lensed maps to 
good accuracy is to synthesise the CMB map with a resolution higher (at least twice larger) 
than the resolution desired for the lensed CMB map and then remap the pixels of the CMB 
map by the deflection angle onto the lens plane. The size of a map is given by the pixel size 
times the number of pixels, with the pixel size being determined by the beam size. 

We then introduced the detector effects, first by convolving the lensed map with the beam 
profile hi = ^2 i exp[—(9f w h mi ) 2 £(£ + l)/(8hi2)/2] and then by adding the map constructed 
from the detector noise power spectrum 1/^Jl/ (Ofwhmi&pixi) }, where the sums are over 
the number of channels i. Finally we deconvolved the resulting map with the beam profile. 

Appendix C: Implementation of the harmonic space estimator 

In this Appendix we describe the implementation of this estimator to reconstruct the con- 
vergence field kq. The optimal estimator is a convolution in harmonic space of the functions 
F\(0) and VeF 2 {d), which are defined as follows 

F X {1) = -=J— " f{£) (CI) 

F 2 (£) = J^—- f(£), (C2) 



C t + N t 



so that 



m = -iAf i J^£-£'F 1 (£-£') (VeF 2 )(£') 
= -iM t l- f d 2 exp[-i£ ■ 0] F 1 (0) V ' F 2 (6) 
= -Af t J d 2 6 exp[-i£ • 0] V [^(0) V e F 2 {6)} . (C3) 

The estimator of the lensing field measures the correlations between the small-scale filtered 
temperature, as encoded in Fi(0), and the temperature gradient, as encoded in VeF 2 {6). 
This result can be intuited by noting that while on small scales the unlensed temperature 
anisotropy is approximately a temperature gradient, the small-scale lensed anisotropy comes 
from disturbing this gradient by the deflection angle from the lensing potential. Note that 
the quadratic estimator is only valid to lowest order in the lensing potential. 
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Figure 8: Geometric relation of the position vectors in harmonic space. 

From the lensed CMB map previously synthesised, we first construct the two scalar fields 
Fi(£) and F 2 (£). We then construct two high-pass filters in real space, one explicitly from 
Fi (£) and the other from the gradient of F 2 {£) 

d 2 £ 



F 1 (0) 
VeF 2 (6) 



(2tt) 2 
d 2 £ 

WT 2 



F^t) e 



it-B 



i£ F 2 {£) e 



m-e 



(C4) 

(C5) 



We can regard the first filter as a weight of the gradient of the temperature map (the second 
filter), resulting the vector field F 3 (0) 



F 3 (0) = F 1 (0)VgF 2 (0). 



(C6) 



We then take the divergence of F 3 (0) and transform into the £ space, defining the scalar 
field F 4 (£) 



F 4 (£) = i£-F 3 (£). 



(C7) 



This quantity is essentially the quadratic estimator of the gravitational potential ip up to 
the total variance Mi 



^(£) = -M F 4 (£). 



(tt 



(See Appendix D for the parametrization adopted to calculate A^.) Finally we compute 

a{£) = i£ $(£) =Af e £[£- F 3 (£)} (C9) 

«o(£) = f 4>(£). (CIO) 



Appendix D: Calculation of the total variance 

In this Appendix we detail the coding of the integrand function for the calculation of 



the total variance Mi given by Eqn. (A13). First we demonstrate the calculation in the 



parametrization (£,£') and then extend the reasoning to the parametrization (£ + ,£^ 
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Since the intervening vectors form a triangle, we adopt for coordinate system an orthog- 
onal grid (i,j) such that the j axis coincides with the direction of £, i.e. £ = £j. The other 
two vectors £' and £" = £ — £' have components along these axes as follows (see Fig. M) 



t = hi + [ h + g ) 3 (Dl) 

£" = -£ 2 i - ( £l - 1 ) j. (D2) 



In this coordinate system we work out the quantities 



02 1)2 



£.£' = i£ 1+ £. £" = -ti x + (D3) 



. 2 

' 2 _ i>2 , / n i l \ c" 2 _ «2 



as well as 

and implement a two-dimensional numerical integration routine over £ x and £2, with £ as an 
input, to compute Mi- 

When we instead work in the coordinates (£+,£_), we note that £ = £ + and £' = (£ + + 
£-)/2 which implies that £" = (£ + — £J)/2. Following the same reasoning as above, we find 
the same component decomposition for £' and £" along the (i,j) axes, only with £ replaced 
by£ + . 
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